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ABSTRACT 

We use three-dimensional direct numerical simulations of the helically forced magnetohydrodynamic equa- 
tions in spherical shell segments in order to study the effects of changes in the geometrical shape and size of 
the domain on the growth and saturation of large-scale magnetic fields. We inject kinetic energy along with ki- 
netic helicity in spherical domains via helical forcing using Chandrasekhar-Kendall functions. We take perfect 
conductor boundary conditions for the magnetic field to ensure that no magnetic helicity escapes the domain 
boundaries. We find dynamo action giving rise to magnetic fields at scales larger than the characteristic scale 
of the forcing. The magnetic energy exceeds the kinetic energy over dissipative time scales, similar to that 
seen earlier in Cartesian simulations in periodic boxes. As we increase the size of the domain in the azimuthal 
direction we find that the nonlinearly saturated magnetic field organizes itself in long-lived cellular structures 
with aspect ratios close to unity. These structures tile the domain along the azimuthal direction, thus resulting 
in very small longitudinally averaged magnetic fields for large domain sizes. The scales of these structures are 
determined by the smallest scales of the domain, which in our simulations is usually the radial scale. We also 
find that increasing the meridional extent of the domains produces little qualitative change, except a marginal 
increase in the large-scale field. We obtain qualitatively similar results in Cartesian domains with similar aspect 
ratios. 

Subject headings: MHD - Turbulence 



L INTRODUCTION 

A fundamental question in solar and stellar physics con- 
cerns the generation of large-scale magnetic fields in con- 
vective spherical shells through dynamo action, which occurs 
on dynamical time scales. A great deal of effort has gone 
into understanding this question by using direct three dimen- 
sional magnetohydrodynamic (MHD) simulations, in Carte- 
sian domains with forced and convective turbulence as well 
as in spherical domains. These studies can be divided into 
four broad groups. The first consists of helically forced tur- 
bulence simulations in Cartesian domains, see e.g., Branden- 
burg (2001); Brandenburg & Dobler (2001). These simula- 
tions in general show large-scale magnetic fields when peri- 
odic or perfect conductor boundary conditions are used, but 
only growing on dissipative time scales, which makes them 
not directly relevant to solar and stellar situations. With more 
realistic open boundary conditions and in presence of shear, 
large-scale magnetic fields are known to develop on dynami- 
cal time scales (Brandenburg 2005). The second group com- 
prises simulations of turbulent convection in Cartesian coordi- 
nates, which have recently shown large-scale magnetic fields 
(Kapyla et al. 2008; Hughes & Proctor 2008). Thirdly, 
forced incompressible turbulence simulations in full spheres, 
mostly relevant to planetary dynamos, have been carried out 
by Mininni & Montgomery (2006); Mininni, Montgomery, 
& Turner (2007). Finally, there is an increasing body of 
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work employing simulations of MHD turbulence in spheri- 
cal rotating shells with convection using the incompressibil- 
ity constraint with either Boussinesq approximations (Oilman 
& Miller 1981; Oilman 1983) or anelastic approximations 
(Oilman & Olatzmaier 1981; Olatzmaier & Oilman 1982; 
Olatzmaier 1984, 1985; Miesch et al. 2000; Brun et al. 
2002, 2004, 2006; Brown et al. 2007). These simulations 
produce mainly small-scale magnetic fields and only insignif- 
icant large-scale magnetic fields with parameters relevant to 
the solar and stellar settings. Relatively stronger large-scale 
(global) magnetic fields have, however, been found in rapidly 
rotating shells (Brown et al. 2007). Also, it has recently been 
shown that in simulations of fully convective stars the energy 
in the longitudinally averaged magnetic field can become lo- 
cally comparable to the kinetic energy (Browning 2008). 

In the present paper we attempt to bridge the gap between 
studies in Cartesian and spherical shell domains by solving 
the MHD equations in wedge-shaped domains of spherical 
shells with helical forcing. In particular we study the effects 
of shape and size of the computational domain on the growth 
and saturation of the large-scale magnetic field. Spherical 
wedge geometries in principle provide an advantage in terms 
of computational resources over both the Cartesian boxes and 
spherical shell geometries usually employed in MHD simu- 
lations in that they strike a reasonable compromise between 
the requirements for spatial resolution and globality. In other 
words, our choice of spherical wedge domains allows in prin- 
ciple higher absolute spatial resolution (i.e., higher number of 
grid points per unit length), thus potentially allowing larger 
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magnetic Reynolds numbers (defined later) to be attained, 
whilst retaining some globality. Alternatively, at a given res- 
olution, we can achieve simulations in a number of wedge 
domains, or in one domain for much longer time, for the cost 
of one simulation in a full spherical shell - this is the approach 
we adopt here. 

In this paper we make a number of assumptions that are mo- 
tivated by the desire to understand the basic concepts of dy- 
namo saturation in spherical geometries instead of providing 
a realistic model of the solar dynamo. Specifically, we con- 
sider here the case of homogeneous turbulence with perfectly 
conducting boundary conditions so as to make contact with 
corresponding earlier work in Cartesian domains. The physi- 
cally more relevant case of open boundary conditions with an 
equator and differential rotation or shear will be postponed to 
future work. 

The paper is organized as follows. In Sect. 2 we briefly 
describe the details of our model and the code used. Sect. 3 
contains our results, where for the sake of clarity, we present 
the results concerning the effects of increasing the domain in 
the azimuthal and meridional directions separately. Sect. 4 
contains our conclusions. Finally, Appendices A and B con- 
tain the details of the helical forcing used, and our extension 
of the Pencil Code^ to non-Cartesian coordinate systems, 
respectively. 

2. THE MODEL 

We solve numerically the magnetohydrodynamic equations 
for the velocity U, the logarithmic density In p, and the vector 
potential A, given by 

DtU=-clV\np+^JxB + F^i,^ + f, (1) 

Dtlnp^-V-U, (2) 

dtA = U X B + r]V'^A, (3) 

where Fvisc = (/^/p)(V^C/ + -|VV-C/) is the viscous force, 
fi is the dynamic viscosity, i? = V x A is the magnetic field, 
J = V X is the current density, po is the vacuum 

permeability, Cg is the velocity of sound in the medium, p is 
the density, ry is the magnetic diffusivity, and Dt = dt + U ■'V 
is the advective derivative. Here f{x, t) is an external random 
helical forcing (the details of which are given in Appendix A), 
satisfying the condition, 

/ • V X / > 0, (4) 

in order to ensure positive helicity injection over the entire 
sphere. Such a model is reminiscent of constant a effect 
spheres that were studied in the early days of mean-field dy- 
namo theory (Krause & Steenbeck 1967). Similar cases rele- 
vant to planetary dynamos have also been studied recently by 
direct numerical simulations (Mininni & Montgomery 2006; 
Mininni, Montgomery, & Turner 2007). 

A sketch of the meridional cross-section of a typical wedge- 
shaped domain used in our simulations is given in Fig. 1 . We 
confine ourselves to simulations in the northern hemisphere. 
However, because there is no rotation, the choice of the co- 
ordinate axis (and hence of the equator) is arbitrary. There- 
fore the physical conditions are the same on either side of 
the equator. The code used for our computations is the Pen- 
cil Code developed by Brandenburg & Dobler (2002) in 

' http://www.nordita.org/software/pencil-code. 



Cartesian coordinates. We have extended the code to allow 
simulations in spherical coordinates. This was facilitated by 
the fact that the Pencil Code was already written in a non- 
conservative form, which allowed the curvilinear coordinates 
to be implemented by replacing all partial derivatives by co- 
variant derivatives, see Appendix B for further details. 




X 

Fig. I. — Schematic representation of the meridional plane of our 
spherical wedge computational domain. We also define 62 to be the 
angle that the other azimuthal boundary makes with the polar axis; 
in this Figure and throughout this paper O2 = 7r/2. 

Guided by the convection zone of the Sun, in the majority 
of our computations the radial extent of our domain is chosen 
to be 0.7 < r < 1.0. We use perfect conductor boundary 
conditions for the magnetic field to ensure that no magnetic 
hehcity escapes the domain boundaries. In Cartesian domains 
(Brandenburg 2001) this is often achieved by assuming peri- 
odic boundary conditions across the boundaries. In our spher- 
ical case this ttanslates to the normal component of the mag- 
netic field B being continuous (and hence zero) across the 
boundary. This implies that the tangential components of the 
magnetic vector potential A must be zero at the boundary. We 
are free to choose the boundary condition for the normal com- 
ponent. Guided by this, we make the following choices at the 
four boundaries of our domain 

dAr 

Ae = A4, = —— =0 (on r = n), (5) 
dr 

A0 = A^ = Ar = Q (on r = r2 and 6» = 6*2 = 7r/2), (6) 
^r = A^ = ^=0 (on^ = ^i). (7) 

There is no particular reason for using = on r = r2 
and not on r = n, and we emphasize that the condition on 
the normal component of A is of no significance for B itself. 
We use stress-free boundary conditions for the velocity at all 
these four boundaries and periodic boundary conditions for 
all the variables along the azimuthal direction. 

3. RESULTS 

Our principle aim in this paper is to study the growth and 
saturation of large-scale magnetic field in spherical wedge do- 
mains. In particular we study the effects of changes in the 
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TABLE 1 

Parameters of the spherical runs. 



3.1. Initial growth phase 

We first summarize our results concerning the growth phase 
_pf the dynamos. We begin with the run with the smallest do- 



Runs 



Grid 



Re Re 



M 



K 



51 32 X 32 X 32 O.Itt O.Itt 

52 32 X 32 X 128 O.Itt n/2 

53 32 X 32 X 256 O.Itt tt 

54 32 X 64 X 32 0.2tt O.IOtt 

55 32 X 256 X 32 85° O.Itt 
86 32 X 64 X 128 0.2Tr Tr/2 
37 32 X 32 X 64 CItt 7r/4 



0.5 
0.4 
0.4 
0.5 
0.4 
0.4 
0.2 



5 14 0.66 

5 12 0.74 

5 12 0.74 

5 12 0.65 

5 12 0.73 

5 11 0.79 

2 4 0.79 



A main size, i.e., SI; see Fig. 2. As can be seen the magnetic 
Q Qg energy starts growing exponentially from t 



0.2r^ and the 



shape and the size of the domain on the resulting large-scale 
fields. For the sake of clarity we do this by studying the effects 
of increasing the domain extent in the 6 and (j> directions in 
turn. We also briefly look at the role of the radial extent of the 
computational domain. Given that the size of our simulation 
domain is different along different directions in different runs, 
we in general have three different length scales, = r2 — ri, 
Lg = r2{02-0i) andL^ = r2 sin 6*2 ((?!)2 -i/)!), corresponding 
to the sizes of the domain in the r, 6 and (p directions respec- 
tively. As an estimate of the characteristic Fourier mode of 
forcing we use k{ = Wrms/f^rms, where M^rms = (W^)^/^ 
is the rms value of the vorticity, W = V x U, and C/rms 
is the rms velocity. Here, angular brackets denote volume 
averages. The characteristic length scale of forcing is de- 
fined to be £{ = 2TT/kf. We then define the fluid Reynolds 
number, magnetic Reynolds number and the turnover time as 

Re = U^ras/vki, Rcm = Urms/Vk[ and T = (C/rmsfcf)"^ 

respectively. Here, u is the kinematic viscosity given by 
ly = n/po where po is the initial density (which is equal 
to the mean density throughout, noting that the mass in the 
volume is conserved). In all our runs r is nearly the same 
and varies from 0.6cs/r2 (run S7) to 0.9cs/r2 (run SI). As 
the dynamo we study is resistively limited (cf. Brandenburg 
2001), the time scale of saturation is the dissipative time scale 



0.10 total magnetic energy reaches the level of the kinetic energy 
0.10 at t » 3r^. This is true of all our runs, since they all start 
^ with the same initial field strength, and they all have the same 
Q growth rate which, in turn, is proportional to U^^^ki, which is 
0.10 also the same for all runs. The growth rate during this expo- 
nential growth phase is given by 



Tr^ = rjkf which is used to non-dimensionalize the time. 

= 7]k\ which we use to normalize the time axes of our 
plots. Here fci is the wavenumber corresponding to the small- 
est length scale in our domain, i.e., ki = 27r/ir for most of 
our runs. The helical nature of the velocity field is character- 
ized by = (W^ • t/)/(Wrmsf^rins)- Wc Start our simula- 
tions with a zero velocity field and a Gaussian random mag- 
netic vector potential such that the amplitude of the magnetic 
field is of the order of 10~^ in units of (/9o/^o)^^^Cs. 

The growth and saturation of the magnetic dynamo is mon- 
itored by the total magnetic energy, Eyi = (B^) / 2/io, and the 
kinetic energy, Ek — {pU'^)/2. We define the large-scale (or 
mean) magnetic field using longitudinal averaging, 



B{r,e,t) 



2nj 



Bd<f>, 



(8) 



over the extent of the domain. The total energy in the large- 
scale magnetic field is then defined by i?Ls(i) = {B'^)/2po. 

A measure of the level of turbulence in our simulations is 
the Reynolds number. Re, given in Table 1, which is about 
5 in all the runs except run S 7 in which case it is 2. For all 
practical purposes there is essentially no inertial range in the 
spectrum of the fluid obtained from our runs. The 
summary of the runs together with their domain sizes, resolu- 
tions and other relevant parameters are given in Table 1 . In the 
following subsections we summarize the results of our simu- 
lations as the domain sizes in the azimuthal and meridional 
directions, and Lg respectively, are changed separately. 



(9) 



which is about 0.1 for all the runs performed here; see Ta- 
ble 1. In all cases decreases (i.e. it is quenched) after 
Em reaches saturation. We note that even after reaching satu- 
ration the field keeps growing somewhat, similar to what has 
been seen earlier in Cartesian domains with periodic boundary 
conditions (Brandenburg 2001), or with perfectly conducting 
boundaries (Brandenburg & Dobler 2002). Both Ek and Hk 
decrease slightly (by less than 10%) after saturation is reached 
for runs S2, S3, and S 6. For other runs the Hk decreases a 
little more (by factors from about 0.8 to 0.6). 




Fig. 2.— Evolution of (U^) (continuous), (B^) (dashed) and (B^) 
(dash-dotted) during early times from run S 1 . Similar exponential 
growth of the magnetic energy is seen in all the other runs. 



3.2. Formation of large-scale magnetic field 

For the smallest domain chosen here, i.e. S 1 which is clos- 
est to a cube {Lr k, Lq k, L^), we obtain results that are 
very similar to those found earlier from Cartesian simulations 
(Brandenburg 2001). The large-scale magnetic field grows, 
reaches a value close to equipartition and then shows a slow 
saturation on dissipative time scales, see Fig. 3. As we are 
using perfectly conducting boundary conditions the growth 
of the large-scale magnetic field is limited by the decay of 
small-scale magnetic helicity. This has been used to model 
the saturation of the magnetic energy of the large-scale field 
(Brandenburg 2001), 



B 



eq 



effcf 



1 _ f.-'^vkl.it-U. 



(10) 
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Fig. 3. — Evolution of kinetic energy (17^) (continuous), magnetic 
energy (B^) (dasiied) and energy in the large-scale magnetic field 

(B ) (dash-dotted) during late times from run SI. The saturated 
value of the energy in the large-scale magnetic field is comparable to 
the kinetic energy. 

Here tgat is the approximate time when the small-scale field 
has saturated, and km is a new effective wavenumber which 
is related to fci and is treated here as a fit parameter that 
is chosen to match the simulation result. We obtain fc.,„ w 
0.7fci. Here B^^ corresponds to the kinetic energy density, so 
B^q/fiQ = (pu^), and is approximately equal to the energy in 
the small-scale magnetic field. Expression (10) fits the data 
from our simulations quite well as shown in Fig. 4. 

The evolution of the large-scale magnetic field follows 
closely the evolution of the total magnetic field after the time 
when the amplitude of the large-scale field has become steady. 
The growth of large-scale structures for this run (run SI) in 
the equatorial plane and the meridional plane are shown in 
Figs. 5 and 6 respectively. Large-scale structures in the con- 
tour plots of magnetic field in the equatorial plane appear as 
early as i « 500 (about 6t,,) and at late times they encompass 
the whole azimuthal extent of the domain. 




B,, t =1000 



B,, t =1990 



0.7 0.6 0.9 1.0 




0.7 D.I 



Fig. 5. — Contour plots of in the equatorial plane of the domain in S 1 at 
different times showing the gradual establishment of a large-scale magnetic 
field. Time is here given inunitsofi?/cs. 
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Fig. 4. — Late saturation behavior of the mean field B for the run SI 
compared with the prediction given by Eq. (10) (dashed). 
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3.3. Effects of increasing the azimuthal extent 



Fig. 6. — Contour plots of in the meridional plane of the domain in SI 
at different time. Time is here given in units of R/cs- 
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To proceed, we begin by increasing the domain size in the 
(f) direction, while keeping the 6 and r dimensions fixed, and 
increase thereby the aspect ratio. The initial growth phase 
remains practically unchanged as we change the extent of our 
domain. We find that the large-scale magnetic field decreases 
as we go to larger domains (by increasing L^), as can be seen 
in Fig. 7. 

To understand the reason for this decrease in the large-scale 
magnetic field we present contour plots of the 9 component 
of magnetic field in the equatorial plane for four different 
domain sizes at later times (see Fig. 8). Notice that as we 
increase our domain size 'cell-like' structures are developed 
along the azimuthal direction, with aspect ratios close to unity. 
Their typical length scale corresponds to and seems to be de- 
termined by the smallest dimension of our domain, which here 
is the radial extent. We checked this by performing a run with 
half the radial extent and found that the characteristic hori- 
zontal scale of the cell structures is decreased accordingly, to 
half the original value. We also find that the length scale of 
these cell structures does not depend on the forcing length 
scale. We verified this by changing the forcing length scale 
along the radial direction and found that this does not change 
these cell structures. The length scale of the cells is also larger 
than the characteristic length scale of the velocity. This is best 
described using Fourier transform along the azimuthal (peri- 
odic) direction, 



-I 



Um{r, 0) = U (r, 0, (p) exp{imq 



27r' 



Br 



i{r,9)= / S(r, 6*, 0)exp(im0) — . 

J ^TT 



(11) 



(12) 



We can define the averaged spectra of these Fourier trans- 
formed quantities as, 

Sl = {\Um{r,6)\')re, = {\B„,{r,e)\')re, (13) 

where the subscript rO denotes meridional averaging. Note 
that S^n = (U^) and J2 = (B^)- We plot in Fig. 9 
both and for the runs SI, S2 and S3, which have az- 
imuthal extents tt/ 10, 7r/2 and TT respectively. We find thatthe 
peak in the spectrum of the magnetic field occurs at the same 
m for runs S2 and S3, showing that the typical characteristic 
scale of the periodic structures does not depend on the </> ex- 
tent of our domain. Note also that the typical forcing scale is 
clearly smaller than this (corresponding to m « 20-40). 

An important question regarding these structures, and hence 
the resulting large-scale magnetic fields, is whether these pe- 
riodic structures are transient and may later merge to form 
structures encompassing the whole domain similar to the run 
SI with the smallest domain size. The characteristic time 
scale over which structures encompassing the whole domain 
form in SI is about 6t^. In Cartesian simulations with mag- 
netically closed boundaries the saturation time is inversely 
proportional to the square of the relevant domain size. Thus, 
by analogy, if the (p extent is doubled the time scale ~ 6t^ 
would become 24t^. Similarly it would take even longer 
for such structures to form in the runs with bigger domain 
size. We have studied a run - run S 7 - in which the azimuthal 
extent of the domain is twice that of SI and have run this 
simulations up to 500t^, without finding any evidence of cells 
merging. This suggests that the periodic structures that we ob- 
serve are at least as long-lived as the duration of our longest 
runs. 



To summarize, our simulations show that the characteristic 
scale of the large-scale magnetic fields found in our simula- 
tions is about the scale of the radial extent of our domain. 
Hence, as we increase an increasing number of periodic 
structures appear along the azimuthal direction, which, in the 
largest domain we have used, is about 10 times larger than the 
radial direction. Therefore the large-scale magnetic field, de- 
fined as a longitudinal average, gives a very small contribution 
in the runs with larger domains. Note that, with this definition, 
the energy of the large-scale magnetic field corresponds to the 
energy in the (axisymmetric) Sq mode. From the plot of the 
spectrum we note that most of the magnetic energy is actu- 
ally concentrated at m = 8, which is the scale of a cell, and 
this mode indeed shows super-equipartition. Also note that 
this mode corresponds to length scales larger than the scale of 
forcing, which corresponds to m « 20 — 40. Hence, instead 
of using the longitudinal average to calculate the energy in the 
large-scale magnetic field we can use the energy in the mode 
m = 8 of the meridionally averaged spectrum of the magnetic 
field. A comparison between these two methods of calcula- 
tion of energy in the large-scale magnetic field is shown in 
Fig. 10. Note that as is increased the large-scale magnetic 
energy measured by the averaged spectrum of the magnetic 
field remains practically constant. 
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Fig. 7. — Normalized energy in the large-scale magnetic field versus time 
for the runs SI, S 2 , S3 and S 7 . As can be seen, as the domain size increases 
in the </>-direction the field decreases (see also Fig. 10). The inset shows the 
same plot but in linear scale for the run S7 which was run more than 10 times 
longer than the other cases. 



3.4. Effects of increasing the meridional extent 

Next we study the effects of increasing the domain size by 
increasing the meridional extent. We find that, as we increase 
the domain size along the 6 direction, the large-scale magnetic 
field shows marginal increase; see Fig. 1 1 . Contour plots of 
the toroidal component of the magnetic field in the merid- 
ional plane for the runs SI, S4 and S5 are shown in Fig. 12. 
Note that, as the domain is increased in the 6 direction, the 
field structure at low latitudes is largely unchanged, while new 
weaker fields are added at high latitudes. However, the high 
latitudes contribute relatively httle to the volume average, so 
the magnetic energy is only marginally increased. As noted 
above, our use of the term high latitudes is defined by our ar- 
bitrary choice of the coordinate axis - see Fig. 1 . However, 
once such a choice is made, the field can only develop subject 
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Fig. 8. — Contour plots showing the typical structure of the magnetic field 
in the equatorial plane as the (f) extent of the domain is increased. From top 
to bottom, plots of the runs S3,S2, S7 and SI. 
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Fig. 9. — Kinetic and magnetic energy spectra 5*^ (dashed line) and 
(continuous line) from runs SI, S 2 , and S 3 . The range of scales 
being forced are shown within the two arrowheads. For clarity the 
spectrum for run SI and S2 are multiplied by a factor of 10~ and 
10^^ respectively. 



10.00 



01 « 

A 

105 

V 



1.00 



0.10 



0.01 




0.0 0.1 0.2 0.3 0.4 



0.5 



Fig. 10. — Two different measures to estimate the energy in the 
large-scale magnetic field plotted against L^, for four different runs 
having different domain sizes along the azimuthal direction. In one 
case the large-scale magnetic energy is estimated by longitudinal av- 
erage (denoted by * in the plot), in the other case (denoted by o in 
the plot) it is estimated by the magnitude of the peak of S^. 
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to the constraints imposed by the geometry of the computa- 
tional domain. 

Again, we have checked that the characteristic scale of the 
cells is determined by the smallest of the three dimensions of 
our domain, by performing a simulation in which the merid- 
ional extent of the domain is half that of run SI. We find 
that the resulting cells again have length scales comparable 
to the smallest scale of the domain which in this case is the 
meridional extent. Furthermore, we have checked that, as we 
increase our domain along the azimuthal extent, the cell-like 
structures in our simulations are independent of the merid- 
ional extent of the domain, provided that the radial scale re- 
mains the smallest. To illustrate this we compare contour 
plots of Bg in the equatorial plane for the runs S 2 and S 6 
in Fig. 13. These two runs have the same azimuthal extent, 
Ltf, = 7r/2, but different meridional extents, viz., Lg — tt/IO 
for the run S2 and Lg = tt/A for the run S 6. As can be seen 
the cell-like structures that appear have the same global fea- 
tures. Finally we have checked that the characteristic length 
scale of cells is independent of the forcing scale by perform- 
ing a simulation in which the characteristic scale of forcing is 
half that of the scale of forcing in run SI. 



Fig. 1 1 . — Normalized energy in the large-scale magnetic field versus time 
for three different runs SI, S 4 , and S 5 . The inset shows the same plot but in 
linear scale. 



TABLE 2 

Summary of the runs, including the extents of the 
computational domains for the cartesian runs. 



Runs 


Grid 


Ly L/z 




Re=ReM 


Hk 


T 


CI 


32 X 32 X 32 


6 6 


0.5 


16 


0.7 


2.6 


C2 


32 X 32 X 128 


6 12 


0.7 


25 


0.6 


3.6 



3.5. Cartesian versus spherical: geometry versus aspect 

ratio 

An important question concerning our results is how to 
differentiate between the effects of geometry (globality) and 
changes in the aspect ratio. To answer this question we need 
to compare our wedge domain simulations with simulations 
in Cartesian boxes with appropriate aspect ratios. To this end 
we perform two simulations in Cartesian coordinates - runs 
CI and C2 - with aspect ratios one and two respectively. Rel- 
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Fig. 12. — Meridional cross-sections of after saturation has been 
reached. From left to right: SI, S4, and S5 




Fig. 13. — Equatorial cross-sections of Bg after saturation has been 
reached from runs S2 and S 6, 

evant parameters for these runs are summarized in Table 2. 
These Cartesian simulations correspond to the runs SI and 
S7 in the spherical wedge domains. The main features of our 
spherical runs are also found in the Cartesian runs. In particu- 
lar, we find that the initial (kinematic) growth rate is the same 
for the runs CI and C2. We also observe formation of cell- 
like structures with unit aspect ratio in the run C2. Figure 14 
gives a summary of our Cartesian simulations showing plots 
of kinetic and magnetic energy versus time. Comparing these 
results with the corresponding plots for our spherical wedge 
runs S 1 and S 7 we observe that the decrease in the large-scale 
magnetic field is similar to those in the Cartesian domains, if 
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Fig. 14. — Normalized energy in the large-scale magnetic field versus time 
for the runs CI, C2. As can be seen as the domain size increases in the 
^-direction, the field decreases in a way very similar to the spherical case. 

the aspect ratios are chosen similarly. Similar behavior has 
earlier been seen in simulations in Cartesian domains with as- 
pect ratios not equal to unity (Brandenburg, Dobler, & Subra- 
manian 2002). What is particularly interesting in our case is 
that, in both Cartesian and spherical coordinates systems, the 
observed cell structures are persistent, with lifetimes larger 
than the duration of our longest simulations which, in turn, 
are longer than the magnetic diffusion time based on (for 
spherical runs) or (for Cartesian runs). For example, we 
have checked that in the case of run S 7 the cell structures re- 
main unchanged for at least as long as 400 dissipative times, 

4. CONCLUSION 

We have made a detailed numerical study of the effects of 
changes in the geometrical shape and size of the spherical 
wedge domains on the growth and saturation of large-scale 
magnetic fields. We have used direct three-dimensional nu- 
merical simulations of helically forced MHD equations using 
random helical forcing. 

For the smallest domain with aspect ratio close to one we 
find dynamo action resulting in magnetic fields on scales 
larger than the characteristic scale of the forcing. The large- 
scale magnetic energy grows to exceed the kinetic energy over 
diffusive time scales, similar to that seen earlier in Cartesian 
simulations in periodic boxes. 

In domains larger in the azimuthal direction the large-scale 
magnetic field organizes itself in cell-like structures in the az- 
imuthal direction. The aspect ratio of the individual cells is 
close to unity. This large-scale pattern in the azimuthal di- 
rection has m = 8. This is determined by the smallest (ra- 
dial) scale in the simulations. [As an aside, somewhat sim- 
ilar behavior was found by Moss, Tuominen & Brandenburg 
(1990) in a study of mean-field dynamos.] To encompass such 
a structure within the computational domain we have to have 
Lij, > (7r/4)r2. Further increases in the size of the domain in 
the azimuthal direction just make the cells repeat themselves, 
tiling the domain along the azimuthal direction, and thus re- 
sulting in very small longitudinally averaged fields for larger 
domain sizes. We note that this implies that extrapolation of 
results from Cartesian box simulations with unit aspect ratio 
to spherical shells can be misleading. We have also studied 
the effects of increasing the size of the domain in the merid- 
ional extent. The resulting large-scale magnetic fields show 
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little qualitative change except for a marginal increase, pro- 
vided Lg > (7r/5)r2; see Fig. 8. Hence the smallest wedge 
shaped domain in which we can expect to observe features 
of simulations in a full sphere must have = (7r/4)r2 and 
Le = (7r/5)r2. 

Furthermore the presence of the cellular structures along 
the azimuthal direction means that the usual employment of 
longitudinal averaging loses much information if used as a 
way to define large-scale magnetic fields. Clearly a possi- 
ble alternative is to define the large-scale magnetic field via 
Fourier transform along the (j> direction of our domain. A 
large-scale magnetic energy defined in this fashion results in 
strong fields of equipartition strength in all the domain sizes 
we have used. We note here that the large-scale magnetic field 
defined using Fourier filtering obeys some of the the Reynolds 
rules only approximately. For example the average of the 
product of an average and a fluctuation vanishes only for infi- 
nite scale separation. This shortcoming may cause some dis- 
crepancies between theory and model, which is however be- 
yond the scope of the present paper. 

In all our simulations we find the cellular structures to be 
long-lived with lifetimes longer than the duration of our sim- 
ulations. This therefore suggests that these structures are not 
transients. In an attempt to compare with mean-field dynamo 
models one must note that the excitation conditions for modes 
with rn > 2 are normally much higher than for m below 2, al- 
though there is a clear trend for this difference to diminish for 
thinner shells (Brandenburg, Tuominen, & Radler 1989). On 
the other hand, anisotropics of the a effect might significantly 
change this. 

It is important to clarify the similarities and differences be- 
tween our results and those of previous studies. Forced turbu- 
lence simulations have been carried out by Mininni & Mont- 
gomery (2006); Mininni, Montgomery, & Turner (2007) who 
also adopted a forcing function in terms of Chandrasekhar- 
Kendall functions - although not random - and used perfectly 



conducting boundary conditions. They included the effects 
of rotation and considered both helical and non-helical forc- 
ings. Their computational domain is a full sphere. They con- 
sidered laminar flow patterns and found large-scale magnetic 
fields to be generated, but the energy contained in the large- 
scale component is generally small compared with the kinetic 
energy. Fully turbulent simulations in spherical shells have 
been studied by Brun et al. (2002, 2004, 2006); Brown et al. 
(2007). These flows are subject to rotation and stratification 
which make them helical. However, the degree of heUcity is 
weak compared to our fully helical forcing functions and a 
broad range of wavenumbers is being driven, so it is difficult 
to identify a well-defined energy-carrying scale. 

Our simulations show the effects of magnetic helicity con- 
servation (see Fig. 4), but the magnetic Reynolds number is 
still rather low, so it may be of interest to repeat such sim- 
ulations at larger magnetic Reynolds numbers. However, it 
is important to run for sufficiently long times to be able to 
obtain full saturation. Obviously, such long saturation times 
are not astrophysically relevant, and earlier work in Carte- 
sian domains gives clear predictions that the constraints from 
magnetic helicity are alleviated in the presence of shear giv- 
ing rise to small-scale magnetic helicity fluxes (Brandenburg 
2005; Kapyla et al. 2008). Allowing for latitudinal differen- 
tial shear motions is therefore one of our next objectives. 
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APPENDIX 

RANDOM HELICAL FORCING IN SPHERICAL COORDINATES 

In this Appendix we briefly describe the helical forcing used in our simulations in spherical wedge domains. We require the 
heUcity of the forcing to be positive at every time-step at every grid point. Furthermore in order to obtain a turbulent state we 
use random forcing which is white-in-time. In Cartesian coordinates this is achieved by using appropriately normalized Beltrami 
waves (Brandenburg 2001); in the spherical case we need to use the Chandrasekhar-Kendall function (Chandrasekhar & Kendall 
1957). Similar forcing functions, although not random, in spherical coordinate systems have also been discussed by Livermore, 
Hughes & Tobias (2007). 

To guarantee positive helicity we demand, foUowing Chandrasekhar & Kendall (1957), 

Vxf = af (Al) 
with a positive a at every point in our computational domain. This in turn impUes that / should have the form 

V X V X / = aV, (A2) 

which, using V • / = 0, becomes 

+ a^/ = 0. (A3) 

Clearly all solutions of this equation are solutions of Eq. (Al) but the converse is not true. To find solutions of (A3) consider a 
scalar function tjj satisfying the Helmholtz equation, 

V'^tp + aV = 0, (A4) 
whose solutions in spherical polar coordinates are obtained in terms of spherical Bessel function and spherical harmonics, 

oo / 

V' = E E ^'(aOi"r(^, </>) exp{iU), (A5) 



Dynamos in spherical geometry 



9 



TABLE 3 

Values of a that satisfy Eq. (A9) used in the 







RIN 


SI. 




m 


I 


ai 






20 


81 


129.011139 


135.938721 


143.325378 


20 


83 


130.880829 


137.703308 


144.992371 


20 


85 


132.771484 


139.489746 


146.681885 


20 


87 


134.682465 


141.297455 


148.393219 


20 


89 


136.613068 


143.125763 


150.125793 


40 


81 


129.011139 


135.938721 


143.325378 


40 


83 


130.880829 


137.703308 


144.992371 


40 


85 


132.771484 


139.489746 


146.681885 


40 


87 


134.682465 


141.297455 


148.393219 


40 


89 


136.613068 


143.125763 


150.125793 


60 


91 


138.562683 


144.974060 


151.879028 


60 


93 


140.530670 


146.841827 


153.652344 


60 


95 


142.516434 


148.728455 


155.445251 


60 


97 


144.519348 


150.633453 


157.257141 


60 


99 


146.538788 


152.556305 


159.087616 


80 


121 


169.644516 


174.748535 


180.321686 


80 


123 


171.805023 


176.847809 


182.341858 


80 


125 


173.971619 


178.958252 


184.375427 


80 


127 


176.143417 


181.079208 


186.421967 


80 


129 


178.319519 


183.209961 


188.481140 


100 


121 


169.644516 


174.748535 


180.321686 


100 


123 


171.805023 


176.847809 


182.341858 


100 


125 


173.971619 


178.958252 


184.375427 


100 


127 


176.143417 


181.079208 


186.421967 


100 


129 


178.319519 


183.209961 


188.481140 






ar) = aiii{ 


ar) + bini{' 


ar). 



where 

(A6) 

Here ji and ni are spherical Bessel functions of the first and second kind respectively and ai and bi are constants determined by 
the boundary conditions. A solution of Eq. (Al) can then be constructed as the sum 

/ = T + S, (A7) 

where 



V X (e^), S 



1 , 



-V X T. 



(A8) 



We wish to confine our forcing to certain bands of length scales and also to randomize it. The characteristic scales of the forcing 
function in the radial, meridional and azimuthal direction are given by a, I and m respectively. As to the choice of boundary 
conditions, we demand that / is zero at the two radial boundaries r = ri and r = r2. The constants a;, 6; and a are then related 

by 

aiji{ari) + bini{ari) = aijiiar^) + hniiar^) = 0. (A9) 

For a particular choice of I this transcendental equation has an infinite number of solutions for a and the ratio ai/bi. A higher 
value of a implies more zeros of the function zi{ar) lies within ri and r2, which in turn implies that the characteristic radial 
scale of zi{ar) becomes smaller. Note that we have periodic boundary conditions along the azimuthal direction, hence the 
non-zero values of m which are possible in our domain depends on the extent of the domain in the azimuthal direction, i.e., 
TOmin = 27r/L0, e.g., mmin = 20 for the run SI. In order to mimic turbulence, we force at the intermediate length scales which 
allows kinetic energy to cascade to smaller scales. Furthermore, we want the forcing to go to zero at the equator. This implies I 
must be odd. The values of a, I and m that we use for the run SI are given in Table 3. We used the GNU scientific library ^ to 
compute the Bessel functions and spherical harmonics in our code. 

THE PENCIL CODE IN SPHERICAL POLAR COORDINATES 

The Pencil Code was originally written in Cartesian coordinates. To use it for our simulations of the compressible MHD 
equations in spherical polar coordinates, it needs to be changed accordingly. In fact, given its modularity, the Pencil Code 
is well suited to be generahzed to any curvihnear coordinate system. We do this by writing the MHD equations in a covariant 
form by replacing partial derivatives by covariant derivatives. We shall illustrate this method by considering the particular case 
of spherical coordinates, which is the one relevant to our simulations here. Let us first consider the divergence of a vector field 
A. In Cartesian coordinates using index notation 

V-A = A„,„, (Bl) 

where a comma denotes partial differentiation. The same operator can be written in any non-Cartesian coordinate system by 
replacing the partial derivative by the covariant derivative denoted by a semicolon thus: 

A«^ = A«,/J-r«^^<, (B2) 



■ http://www.gnu.org/software/gsl/. 
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where F^^ are the connection coefficients obtained from the metric corresponding to the coordinates chosen. In the case of 
spherical coordinates the metric takes the form 

/I \ 

9a0={o r-' . (B3) 

We shall write the covariant derivatives in the non-coordinate bases, by defining a new triplet of coordinate differentials 

df = dr, d^ = r d^, and d4> = r sin 6 d(p. (B4) 
In these bases the connection coefficients take a particularly simple form, 

r%^ = r^^ = -r%^ = -r'<ja-iM (B5) 

r^^ = -r%^=cot0/r, (B6) 

with all other cormection coefficients being zero. This simplification makes this non-coordinate basis particularly appealing for 
numerical simulations. For example, for the divergence of a vector A we obtain 

V ■ A = Aa-a = ^a,a + 2r~'^ Af + r"^cotMg-, (B7) 

where 

A&,& = drA^ + -dgAg + -^d(f>A.. (B8) 
r " rsmO ^ 

Note that in the non-coordinate basis the metric tensor is the Kronecker delta and so the covariant and contravariant components 
of a tensor are one and the same, hence in the above expression we have not distinguished between them. As in Eq. (B7) 
any vector differential operator in curvilinear coordinate system can be written as the sum of two parts: the first involving the 
vector operator in the Cartesian form with added scaling factors r^^ and (r sin 9)~^, and the other part involving the connection 
coefficients. The modular feature of the PENCIL CODE then plays an important role since the derivatives in PENCIL Code 
are computed in a separate module, and all we need to do to adapt the Pencil Code to any non-Cartesian coordinate system 
is to change this derivative module by adding the scaling factors corresponding to the coordinate system chosen. The vector 
operators, e.g., divergence, curl, Laplacian etc, are then calculated in a different module which uses the derivative module. The 
parts which depend on the connection coefficients are added to this module. The other minor changes to the code involves coding 
new boundary conditions and new modules to calculate volume averages. All these changes are now part of the public release of 
the code. 

For completeness, we list here the expressions for the most conmionly used vector differential operators in our code. For the 
curl of a vector field A we have 

V X A = I A,.^ - A^., I = I A,j^ - A^ - I + I -r-K\;, \ . (B9) 

For the advective operator we obtain 

(u • VA)r = UfAf^f + UgA. g + - u^A^ - r"^ u^A^ 

(u • VA)^ = UfAj^ ^ + UgA^j^ + u^Ai,^^ + u^Af - r-^cotOu^A^ 

(u • VA)^ = u^A^ f + ugA^ g + u^A^^^ + u^Ar + r-^cot^ u^A^. (BIO) 
To calculate the second order differential operators we need the expression for second order covariant derivative given by 

A ~ — A ■ — r*--4 r*- A - - 

~ &P ~ ~ r'^Q,^ A- J + r'^a/j,r Ac, — r'^^- ^a,* + ^o- (Bii) 

For example the Laplacian of a scalar field ^ is given by 

2 , cot6' 
r 

and the grad div operator takes the form 





A* = E^.^ = (a^*)^ + + — (B12) 



^^d,df + 2r~^ Af^f + cot6Ag ~ - 2r~'^ Af - cot9Ag^ 
VV ■A={ + 2r-i A- g + r'^ coiQAg g - r"^ ^Y^-'^QAg \ . (B13) 
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